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ABSTRACT 

This report presents a detailed operating manual for a grid-generation program that produces three- 
dimensional meshes for advanced turboprops (ATP). The code uses both algebraic and elliptic partial 
differential equation (PDE) methods to generate single-rotation and counterrotation II - or C - type 
meshes for the z - r planes and H - type for the z — $ planes. The code allows easy specification 
of geometrical constraints (such a s blade angle, location of bounding surfaces, etc.) and mesh-control 
parameters (point distribution near blades and nacelle, number of grid points desired, etc.), and it has 
good runtime diagnostics. This report provides an overview of the mesh-generation procedure, a sample 
input dataset with detailed explanation of all input, and example meshes. 

INTRODUCTION 

Three-dimensional computational fluid dynamic codes require grids with suitable resolution and smooth- 
ness, These qualities are especially difficult to maintain because of the taper, twist, and sweep of realistic 
turboprop geometries. Grid resolution near shocks and in shear layers allows for complex fluid physics to 
be captured. Smoothness of the grid prevents truncation errors, arising from the discretization of the flow 
equations, from dominating the solution. 

This report describes a code for generating three-dimensional grids for advanced turboprops. The code 
uses mainly algebraic methods such as (1) Ferguson’s parametric cubic curves, (2) cubic splines, and (3) cubic 
curves to generate meshes. Elliptic partial differential equations (PDE’s) are used for further smoothing in 
certain cases. The resulting meshes, regardless of whether II - or C-type, have in common a two-dimensional 
axisymmetric mesh (common z and r coordinates) , only varying in the 0 coordinate within the passage. 
The resulting meshes may be written in a variety of formats to conform to the input required for different 
flow solvers. For example, for the Adamczyk average-passage code (ref. 1), a counterrotation mesh can be 
written so that eacli blade row has its own H-mesh which describes a full passage from inlet to exit with 
each respective mesh accounting for its own blade thickness while applying zero thickness to the neighboring 
blade row. 


GEOMETRY AND MESH-GENERATION DESCRIPTIONS 

The generation of computational meshes for advanced turboprops consists of a number of steps. First the 
user-supplied nacelle and blade shapes, as well as geometrical and mesh quality constraints are read. From 
this input, the common two-dimensional axisymmetric mesh is created by using a series of parametric cubic 
curves, cubic splines, and, in the case of C-type meshes, by limited use of PDF’s for further smoothing. 
Second, the resulting two-dimensional meshes are smoothed by using equal-weighted averaging. Third, the 
0 - coordinate surfaces are generated, based on the 0 coordinates of the given blade and the user-specified 
blade angle. Finally, the mesh is interpolated to the periodic angle and written subject to user requirements. 

The geometry for the blade and spinner is provided by the user. The spinner/hub shape is input as ( z } r ) 
pairs with z increasing monotonically from the spinner stagnation point or inlet surface, downstream to the 
outflow boundary. The data is in the same scale as the blade geometry. 
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The blade definition is given as blade sections specified on planes perpendicular to a stacking axis which 
is also the pitch-change axis (PCA). The sections (stations) are specified from near the hub to the tip with 
the distance along the stacking axis (rb) increasing monotonically. At a given station, the section is located 
with the local twist angle A/? (dbet) based on the blade angle at 75-percent blade radius /? 3 j A (beta34), 
leading-edge alignment (lea), face alignment (fa), and chord length (cl). The section shape is given in a 
normalized coordinate system (x/cl, ys/cl, yp/cl) with points on the suction and pressure surfaces from the 
leading edge to the trailing edge. Figures 1 and 2 describe this geometry. 

The preceding parameters aTe used to convert the blade sections to a cylindrical coordinate system. Then, 
user-specified point distribution (with respect to uniform unit spacing) and the following hyperbolic tangent 
distribution are used to generate a distribution of normalized points . Let 
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where / is the number of points and As*, As 2 are the fractions of uniform spacing. 
Solving the following nonlinear equation for <5 


gives the unit-segment point distribution as 
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where 




and £ varies from 0 to I. See reference 2. 


By using this normalized distribution, the points on each blade section are redistributed axially. The 
sections are then converted to constant radial sections. The initial curve designs for the two-dimensional 
mesh are done with Ferguson’s parametric cubic curves (ref. 3). These curves are defined as follows: 


r = v(n) — uq T + ^ 3a 3 

where 0 < u < 1. 

To determine ag, a^, a 2 > and a 3 , specify the values of r and dr/ da at both ends of the segments. Now, the 
a’s in terms of the r’s may be obtained in matrix form: 
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The derivatives i*(0) and f(l) are proportional to the unit tangent vectors T(0) and T(l) at the ends of the 
of the line segment under consideration. We may write 


r(0) = a o T(0) ; i-(l) = a,T(l) 

By varying the tangent vector lengths, different-shaped curves are easily generated (fig. 3). These curves are 
used to smoothly splice the z~r curves from the blades to each external boundary. Now the point distribution 
on the blades and the hyperbolic tangent distribution function generate a smooth distribution of points away 
from all solid surfaces. Cubic splines along with the new unit point distributions are used to obtain a very 
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smooth variation of z and r coordinates everywhere. In order to further smooth the z and r coordinates, an 
equal-weighted smoothing is applied. In C-type meshes, the redistribution scheme may cause a high degree 
of skewness in the cone region (the rounded region up to the leading edge of the blade). Elliptic PDE’s are 
used to smooth this area (ref. 2). The basic differential equations under consideration are 

V 2 f = P 

and 

V 2 r; = Q 

where £, rj are the coordinates in the uniform computational domain, and P, Q are called the “control 
functions’ 5 . These functions allow the user to control the placement of the grid points in the desired region. 
However, in this report, the equations solved were 


V 2 £ = 0 


( 1 ) 


and 



( 2 ) 


Hence, the point distribution is controlled only in the 7? direction. However, the solution of equations 1 and 
2 will yield a distribution of £ and ?/ which is known to be a uniform rectangular two-dimensional plane 
(fig. 4). Therefore, inverting equations 1 and 2 in order to change the role of the dependent and independent 
variables gives 


2#i2r -{- i 7 1 1 ^77 J7 — -guQi'ti ( 3 ) 

where r is a vector of the physical coordinates z, r and 

<7 li = V + r f 2 


<712 = ^z J} + r^rjj 


2 * 2 
9 22 = Zrj + TV 

and 


G = (H r r) - Wt) 2 

is known as the square of the Jacobian of transformation. The forcing function may be arbitrary, and the 
one chosen was 

q _ ~(2 + (q - 1) * log(x)) * log(x) 

1- + {l - 1) *log(x) 

This function will allow the user to increase the grid clustering near the 77 = 0 boundary with increasing x 
^ 1) • Equations 3 are first discretized on the computational plane by central differences and then solved 
with a block relaxation scheme where the block was a line in this case (Line SOR). 

Two appendixes are provided to guide the user in the proper use of this program. Appendix A lists and 
explains, in detail, all inputs required; appendix B provides a typical input dataset. 


RESULTS 

Figures 5 to 8 show different views of a counterrotation H-mesh generated by using the attached example 
dataset with JNOSE changed to 8. Figure 5 shows a view of the z — r plane; notice the effect of the clustering 
factors near the nacelle, tip, and at the edges of the blades. Figures 6(a) and 6(b) show the effect of the 
clustering factors (DZLE, DZTE, DZLE2, and DZTE2) on the blade point distribution. DZLE and DZTE 
were changed to 0.5 and 1.0, respectively, for the front blade. Their efFect can be seen clearly in figure 6(b). 
Figure 7 shows a portion of a full passage via a radial cut at the hub, with the axial index starting just a 
few gridlines ahead of the front blade and finishing a few gridlines aft of the rear blade. 
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Figure 8(a) shows the full passage view for a single blade case via a radial cut at the hub. In this case, 
DTITBL was 0.2; DTHBL was then changed to 1.0. The effect of this change is seen in figure 8(b). Figures 
9(a) and 9(b) show a constant axial cut at the inlet and at approximately midchord of the blade, respectively. 

Figure 10 shows an example C-mesh. To generate this mesh, the nacelle input was changed to have no 
sting, CONOPT was changed to True, JCONE was set, and CK was set to 1.1. 

Finally, Figure 11 shows the surface meshes for a swirl-recovery vane (SRV) grid, a cruise missile (CM) 
and an advanced single-rotation turboprop (ATP). 

CONCLUSIONS 

An algebraic and elliptic partial differential equation mesh generator has been developed for generation of 
three-dimensional meshes for advanced turboprops and other similar geometries. H— and C— type meshes may 
be generated for both single-rotation and counterrotation cases. The program has good runtime diagnostics, 
and the user may easily modify the mesh-control parameters. 
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APPENDIX A 

USER’S MANUAL FOR THE ADVANCED TURBOPROP MESH GENERATOR 

The mesh generator is coded in FORTRAN 77 with 34 subroutines and 2 external functions. The code 
has been run on a wide variety of computer systems with no coding changes. Figure 12 shows a flow chart 
of the main program structure. 

The nacelle and blade geometries, the mesh size and spacing parameters, and the controls for miscel- 
laneous options are all read from FORTRAN UNIT 7 in the following formats: (1) The majority of real 
number fields are read in F10.5 format, while a few are read in free-format (F-F), (2) integers are read using 
an 15 format, and (3) logical fields are in read L5 format. In order for the program to operate, the input file 
must be named “fort. 7”. The uppercase lines shown in the following description of the sample input serve 
as comments and are skipped over by the program; however, they must be included in the input for proper 
operation. All actual input lines are listed in their lowercase alphanumeric form. An explanation of all input 
follows the sample. Appendix B shows a sample input dataset for a counterrotation case. 


title There may be as many lines of text in the title section as the user desires. The only restriction is 
that the following line with “ PT10” must exist. Note that the leading space in “ PT10” is significant. 
The title is read with an “A” format; the maximum length of each line is 80 characters. 

PTIO PT20 PT30 PT40 PT50 - (5L5) 

ptlO pt20 pt30 pt40 pt50 These are logical flags which are read with a (5L5) format. The print flags 
(following) print all the coordinates (r, r, 0) ; ifPT30or PT40 has been invoked, ignore the 0 coordinates. 
Until subroutine EXTEND has been called, only the 6 on the blades are correct. When any flag is 
Time, the program will print the following: 

ptlO Print all input, as well as the blade geometry or geometries (if two blade rows) in normalized 
(diameter =z 1.) cylindrical form; see subroutine MEREAD. 

pt20 Print blade(s) after correct axial distribution is applied and the stations are at constant ra- 
dius; see subroutine ADJBLD. Also print the fictitious base coordinates; see subroutine FICBAS. 

pt30 Print the regions comprising the rest of the surface through the blade(s); see subroutines 
FORBLD and MIDBLD (used only when generating two blade rows), OUTREG, and CONREG. 

pt40 Print the blade surface after radial distribution is applied; see subroutine RADJST and then 
after the surface not on the blade(s) is smoothed, see subroutines AVGSRF and AVG. 

pt50 Print the entire, final mesh after all other processing; see subroutine MWRITE. This pro- 
duces hundreds of pages of mostly useless output. Use this option only when necessary. 

COUNTER-ROTATION - (L5) 

cropt Logical counterrotation option flag. Use an L5 format to read. If Ti-ue, then two blade geometries 
will be read from the input, first the front and then the rear. If False, then single rotation is selected . 

NACELLE GEOMETRY 

NUMBER OF INPUT STATIONS FOR NACELLE - (15) 
nn Number of points specified on the nacelle. 

zn(l),rn(l) to zn(nn),rn(nn) The actual nacelle geometry, from the nose or any extension, such as a 
sting for an II-mesh (1) back to somewhere behind the blade (nn). The pitch-change axis (PCA) of 
the (front) blade is used as the axial z origin; upstream of the PCA (toward the nose) is negative. The 
rotation axis is, of course, the origin radially r. The nacelle geometry should be in the same scale as 
the blade geometry(s). The nacelle geometry is read in a 2F10.5 format. 
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BLADE GEOMETRY 

BL.DIA. BETA34 NR NC - (2F10.5.2I5) 


diam The diameter of the blade. Everything is scaled by the program so that the working diameter of 
the (front) blade is 1. in the output mesh. 

beta34 The value of beta at 75-percent radius. Beta is the angle, in a given plane normal to the PCA, 
between the chord line and the plane of rotation. 

nr Number of radial stations in the input. This can be different for the two blade rows, if the user is 
generating a counterrotation grid. The input stations should not include a fictitious station of data at 
rb = 0 (see following for definition of rb). Station nr is at the tip. 

nc Number of stations along the chord line at each blade cross-section in the input. Station 1 is the 
leading edge and station nc is the trailing edge. 

RB DBETA LEA FA CL - (5F10.5) 

X YS YP X YS YP - (10X,6F10.5) 

rb The distance from the rotational axis to the plane containing the current blade section. This should 
be in the same units as diameter. 

dbet The angle, in degrees, to be added to beta34 to get the value of beta for the current station, see 
bcta34. 

lea Leading-edge alignment. The distance along the chord line from the point on the chord line closest 
to the PCA to the leading edge; lea is positive if the point closest to the PCA is aft of the leading 
edge. Same units as diameter. 

fa Face alignment. The perpendicular distance from the PCA to the chord line of the current station; 
fa is positive if the PCA is on the face (pressure) side. Same units as diameter. 

cl The chord length (from leading to trailing edge) at the current station. The following blade coordi- 
nates are nondimensionalized by this value. Same units as diameter. 

x(i),ys(i),yp(i) The blade shape at the current station. Where x is the fractional chord wise distance 
(0. (Pleading edge to 1.0:trailing edge) and (ys, yp) are the distances perpendicular to the chord line, 
normalized by the chord length. The blade section may be input with a finite thickness at the leading 
and trailing edges; the program will, however, force the thickness to be zero at the leading and trailing 
edges. The (face/lower/pressure) surface is defined by yp, while the (camber/upper/suction) surface 
is defined by ys. Note that data for two x values are on each input line. 

MESH PARAMETERS 

ZUP ZDN RINF - (3F10.5) 

zup The upstream boundary of the mesh, measured axially from the PCA in fractions of a diameter. 
It is in a scale where the blade diameter is 1.0; hence, if the input blade diameter changes this should 
stay the same to get the same relative location. It has to be negative. Note that for an IT-mesh, zup 
must not be beyond the upstream limit of the nacelle input, 

zdn The downstream boundary of the mesh. Scaled the same as zup. If it is beyond the downstream 
limit of the nacelle input, the nacelle will be extended straight back (at the same radius) from its last 
station. 

rinf The radial boundary of the mesh. If there is no outer region, .i.e., kmax — ktip (see kmax), this 
will be ignored, and 0.5 (the normalized tip radius) will be used instead. It is scaled the same as zup 
and has to be > 0.5 when used. 
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CONE FLAG - (L5) 


conopt This is the cone option flag. If True, then the cone region (if present; see jcone) will have a 
curved outer boundary (C-mesh nose). If False, then the region will have a rectangular boundary 
(II -mesh). 

JCONE JLE JSL JTE JMAX KTIP KMAX LMAX - (815) 

jcone The location of cone region. J is the index in the axial direction; jcone specifies the number of 
axial stations in the cone region. If it is zero, there will be no cone region; otherwise, it should be at 
least 3. 

jle Index of the leading edge of the blade. The region from jcone to jle is called the “fore region”, so jle 
must be greater than jcone. 

jsl Index of the shock line clustering on the (front) blade. If this is zero, there will be no clustering. See 
dzsl and cfsl. 

jte Index of the trailing edge of the (front) blade. If generating a single-rotation grid, the region from 
here to jmax is the aft region; otherwise, from here to jle2 (leading-edge index of (second) blade row) 
is the midregion. 

jmax Index of the downstream limit of the mesh (z = zdn). J is the axial index. 

ktip Index of the tip of the blade. K is the radial index. The region from 1 to ktip is the region from 
the rotation axis to the blade tip. 

kmax Index of the outer boundary of the mesh. From ktip to kmax is the outer region (excluding the 
part upstream of jcone, which is in the cone region). If kmax = ktip, there is no outer region and the 
tip radius will be the maximum radius instead of rinf. 

Imax Number of stations in the circumferential direction. 1 increases as 0 increases. 

NBLADE - (15) 

nb The number of blades on the front row. 

LOCATE ZLEINF ZTEINF - (L5,2F10.5) 

locflg Logical flag which, when True, will use zleinf and zteinf for the specification of the leading- and 
trailing- edge extensions at the outer boundary. 

zleinf If locflg is True, this specifies, in the units used for zup, the axial location of the extension of the 
leading edge at rinf. It must be less than zdn and greater than zup. If there is a cone region, it also 
should be greater than the scaled-down location of the nacelle nose (zn(l)/diam). 

zteinf Same as zleinf, but for the trailing edge. Must be greater than zleinf and less than (zn(nn)/diam). 
DZLE DZSL DZTE CFSL DRNAC DRTIP DTHBL - (F-F) 

dzle The mesh-clustering factor in the axial direction at the leading edge of the (front) blade. This is 
the fraction of an even spacing on the blade. If it is 1.0, the leading edge will have no clustering, while 
if it is 0.1, the mesh will be approximately 10 times denser at the leading edge (compared to uniform 
unit spacing). 

dzsl Clustering factor for the shockline on the (front) blade, if jsl is not zero. See cfsl. 

dzte Clustering factor for the trailing edge of the (front) blade. 
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cfsl If the shock-line clustering is selected (jsl is not zero), the location of the clustering on the blade 
is given, cfsl stands for chord fraction at the shock line; consequently, if the chord fraction is 0.2, the 
clustering will be one-fifth the distance from the leading edge back to the trailing edge (closer to the 
leading edge). 

drnac Clustering factor in the radial direction at the nacelle surface. In the cone region this also gives 
a compression axially in front of the nose. 

drtip Clustering factor in the radial direction at the blade tip. It is phased out in the fore region so 
that it has no effect in the cone region. 

dthbl Clustering factor in the circumferential direction at the blade surfaces (and their axial and radial 
extensions). 

ZREAR REFL - (F10.5,L5) 

zrear If the user is generating a counterrotation mesh, this specifies the location of the rear PCA with 
respect to the front. Same units as diameter. This has to be greater than 0. 

reflag Logical flag for counterrotation. If Time, the rear blade geometry will be reflected about 0 == 
0. before being used. If the front and rear blades have the same shape (i.e., identical input), only 
reflected, the front geometry can be copied and this flag set. If the ys (suction side) input values are 
greater than the yp values, then this option should be used. Plots of the resulting mesh can also be 
used as guides to determine when to use this option. 

The rest of the input is very similar to the input for the front blade. Some of the variables have a “2” ap- 
pended to their names which indicates that the variables pertain to the rear blade row. The only difference 
that should be noted is that the rear blade is scaled down by the front blade diameter for consistency of 
units. This means that the rear blade tip will not necessarily be at r/d — 0.5. The user is cautioned that at 
this time the option to specify the leading- and trailing- edge extension locations for the rear blade should 
not be used ( it may cause problems); therefore, set locfl2 = False. From the input it may be interpreted that 
the number of blades on the rear row and the number of theta stations in a passage may be different than the 
number from the front. Please note that this is not a general change, because the subroutine MWRITE is 
now very much flow-solver specific. This change is valid only for the Adamczyk flow-solver format because 
it will output two mesh files; one with thickness for the front row and no rear blade thickness and the other 
with no thickness for the front row and with thickness for the rear blade. Sec the flag IPRNT, 

JNOSE - (F-F) 

jnose If an Il-mesh is being generated, then an axial index (jnose) value for the location of the stagnation 
point is read. This helps maintain the original nacelle shape as specified by the input, if jnose is specified 
as 0, then a value of jnose is computed as half of the leading-edge index (jle). If jnose is greater than 
0, then it is taken to be the index desired ( must be greater than 3). If jnose is specified as less than 
0, then the points are distributed axially from j=l to jle, and no stagnation point is preserved. This 
is important for certain types of nacelle bodies that do not have a stagnation point as such. One can 
also set jnose less than 0 for bodies that do have a stagnation point, however, the stagnation point will 
be faired out. if jnose is set greater than or equal to 0 for a body with no stagnation point, the code 
will not run properly. 

CK C-MESII POISSON SMOOTHING PARAMETER ; < 0 NO SMOOTHING) - (F-F) 

ck Used only for a C^type mesh. This is a factor used to pack points above the tip. Values of ck should 
range from I. to 2. Since the Poisson smoothing only operates from ktip+1 to kmax and j — 1 to 
jeone-f 1, the effect of higher ck is to pack more points near ktip+1. A value of ck = 1 will give a 
solution of the Laplace equation. If the mesh appears to have skewness near the leading-edge index 
above the tip without smoothing, then turn on the smoothing; try an initial value of ck = 1.15. Higher 
values will cause more line contraction near the ktip-f I line. 
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IPRNT MESH OUTPUT WRITE FLAG - (F-F) 


iprnt - This flag controls the format of the mesh as it is written. IPRNT — 1 will write a mesh formatted 
for NASPROP-E (ref. 4), IPRNT = 2 writes a mesh formatted for the DENTON code (ref. 5), and 
IPRNT = 3 writes a mesh formatted for the ADAMCZYK codes (ref. I). See subroutine MWRITE 
for actual write statements. 

RAT (CIRCUMFERENTIAL PACKING RATIO; USED IF > 1) - (F-F) 

rat This a circumferential clustering factor. If rat is greater than or equal to 1 then rat will be used for 
the clustering in the 0 direction. This will overwrite any efTect of the dthbl factor given previously. A 
value of rat equal to 1 will give uniform spacing; even small increases in rat (i.e. 0.1 or less) will change 
the packing considerably. In order to use rat, the number of desired mesh points in the circumferential 
direction must be odd; i.e., lmax has to be odd. However, if Imax has been specified as even, the 
program will change lmax to odd. A message informing the user of the change will be written to 
UNIT 6. See dthbl which also does packing in the circumferential direction. When dthbl is specified, 
a hyperbolic tangent stretching is performed; using rat gives a stretching based on series summation. 

OUTPUT 

UNIT 6 is the debugging output; the print flags (ptlO, pt20 , etc.) control what is printed, and informational 
messages (whether error or otherwise) are written here. Also written to UNIT 6 is the header page and any 
text that the user placed in the title portion of the input file. The output should be self-explanatory. UNIT 
11 (and perhaps also UNIT 12, if generating a counterrotation mesh for the Adamcyzk flow-solver) is the 
final output (mesh). Please note that the final output file or files are always named ‘Tort. 11” and ‘Tort. 12”. 
Files by these names should not exist before the mesh generation process because the program will abort 
prematurely. The output files are unformatted. The first record contains the size and location parameters 
and has two forms depending upon whether it is a single-rotation or counterrotation mesh. The actual mesh 
coordinates are then written. See subroutine MWRITE for the actual write statements since they will vary 
depending upon the flag IPRNT. However, all write so that the z, then r, and finally, the 0 (in radians) 
coordinates are written. 
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APPENDIX 13 

SAMPLE INPUT DATASET 
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Figure 5. * z — r plane of counterrotation II -mesh. 



(a) DZTE = 0.5 


(b) DZTE = 1.0 


Figure 6 - Effect of clustering factors on blade point distribution. 
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Figure 7. - Constant radial cut at hub (counterrotation case). 



(a) DTIIBL = 0.1 (b) DTIIBL = 1.0 

Figure 8. - Effect of DTIIBL on circumferential point distribution. 
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At inlet. 


(h) At midchord. 


Figure 9. - Constant axial cut. 



Figure 10. - z — r plane of typical C-mesh. 



(b) Cruise missile (CM) 


Figure 11. - Surface grids. 
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(c) Advanced turboprop (ATP) 


Figure 11. - Concluded. 
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